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Abstract 

Using spin-dynamics techniques we have performed large-scale computer sim- 
ulations of the dynamic behavior of the classical three component Xy-model 
(i.e. the anisotropic limit of an easy-plane Heisenberg ferromagnet), on square 
lattices of size up to 192^, for several temperatures below, at, and above 
Tkt- The temporal evolution of spin configurations was determined numeri- 
cally from coupled equations of motion for individual spins by a fourth order 
predictor-corrector method, with initial spin configurations generated by a hy- 
brid Monte Carlo algorithm. The neutron scattering function S{q, u) was cal- 
culated from the resultant space-time displaced spin-spin correlation function. 
Pronounced spin-wave peaks were found both in the in-plane and the out-of- 
plane scattering function over a wide range of temperatures. The in-plane 
scattering function S^^ also has a large number of clear but weak additional 
peaks, which we interpret to come from two-spin-wave scattering. In addition, 
we observed a small central peak in S^^, even at temperatures well below the 
phase transition. We used dynamic finite size scaling theory to extract the 
dynamic critical exponent z. We find z = 1.00(4) for all T < Tkt, in excellent 
agreement with theoretical predictions, although the shape of S{q, u) is not 
well described by current theory. 



1 Introduction 



The two-dimensional XF-model is one of the 'special' models of magnetism. It undergoes 
an unusual phase transition to a state with bound, topological excitations (vortex pairs) 
but no long range order [0]. The (three component) XY-model may be viewed as a 
special case of the anisotropic Heisenberg model in which the coupling between the z- 
components of spins vanishes. It has static properties which are similar to those of the 
"plane rotator" model, in which the spins have only two components. The static properties 
which of both models have been determined by numerical simulation p|, ^, ^ and found 
to be consistent with the predictions of the Kosterlitz-Thouless theory. For example, 
the susceptibility shows an essential singularity instead of a power law divergence, and 
computer simulations show that vortex pairs unbind at Tkt- The model is critical, i.e. it 
has infinite correlation length ^, at all temperatures T < Tkt- The spin-spin correlation 
function decays algebraically with distance for all T < Tkt, but with a power rj which 
varies with temperature. 

The dynamic behavior of the model should be governed by the dynamic critical expo- 
nent z, which describes the divergent behavior of the relevant time scale 0, i.e. r oc 
for all T < Tkt- Recently, finite size scaling for critical dynamics in the neutron scat- 
tering function has been developed [Q, and successfully applied to the study of a 3-dim. 
Heisenberg ferromagnet. The XF-model has true dynamics which can be determined by 
integrating the equations of motion for each spin; its critical dynamics has been studied 
theoretically [|^ ||, |^, ^ with different predictions for the nature of the neutron scattering 
function. In contrast, the "plane rotator" model does not possess equations of motion; it 
has only stochastic, i.e. purely relaxational time dependence, which has also been exam- 
ined by Monte Carlo simulation [Q. Note, however, that different dynamic exponents are 
expected for stochastic and for true dynamics 

In this paper we present the first large scale simulation study of the true dynamic 
behavior of the XF-model. Great care was taken to ensure that statistical as well as 
systematic errors were both well understood and quite small. An earlier, much less com- 
plete study |Tl| indicated a rich structure in the neutron scattering function which was 
not adequately described by theory. Our model is defined by the Hamiltonian 

n = + sfs]) , (1) 

nn 

where Sj is a three- component classical vector of length unity and the sum is over all 



nearest neighbor pairs. The equation of motion of each spin is [11 



— Si = Sj X Heff, (2) 

where 

Heff = -JJ2 {Sj^o: + S]ey) , (3) 

nn 

and e^^ and By are unit vectors in the x- and y-directions respectively. Eq. (|^) represents a 
set of coupled equations and can be integrated numerically. (The plane rotator model has 
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a hamiltonian of the same form as eq. (|I]), but since the vectors have only two components, 
an equation of motion cannot be defined in the same way as for the XF-model.) 

In section 2 we define the neutron scattering function, provide dynamic finite size 
scahng equations, and summarize analytical results. Section 3 describes the details of 
our simulations. We present our data and subsequent analysis in section 4, and draw 
conclusions in section 5. 

The results of our dynamic study prompted us to reexamine the static properties of 
the model. In order to obtain more reliable values for the critical temperature and the 
static critical exponent r/, we carried out static Monte Carlo simulations; our results are 
presented in the appendix. 



2 Background 

2.1 Neutron Scattering Function S{q,u) 

The neutron scattering function S{q, u) (also called the dynamic structure factor) is an 
experimental observable and is fundamental to the study of spin dynamics. It is defined 
^ for momentum transfer q and frequency transfer u as the space-time fourier transform 

r+oo rlf 

S^\ci,uj) = Ve^'i-^'-'-) / e''^'C''\r-r',t)— (4) 

of the space-displaced, time-displaced spin-spin correlation function 

C"=^(r-r',t-t') = {s'At)SUt')) , (5) 

where k = x,y, or z is the spin component, displacement r is in units of lattice spacings, 
and the angle brackets (...) denote the thermal ensemble average. Note that in the 2-dim. 
XF-model, (^S^{t)'^ = for all components k = x,y,z. The equations of motion (|^) are 

time reversal invariant, therefore C'^'^(r — r',t — t') is symmetric in t and t', and S{q, u) is 
real- valued. 

The neutron scattering function generally depends on the correlation length ^ and 



may be written in the form |12 



-ffcfc/„ , A ^ ckkf„\ fkk I ^ 



where 5'|'^(q) is the total intensity given by 

/oo 
Sl\q,u;)du, (7) 
-oo 

and /^'^ is a normalized shape function, f'''^{x,q^,^)dx = 1. The characteristic fre- 
quency u!^{q,$,) is a median frequency determined by the constraint 



= / ^Sl\q,uj)duj. (8) 
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In dynamic scaling theory it is assumed that the median frequency uj^{q,C,) is a homo- 
geneous function of q and ^, i.e., 

c^™(q,0 = (9) 

where z^'' is the dynamic critical exponent and fl^'' is another shape function, and that the 
function f^'^ depends only on the product of but not on q and ^ separately. Therefore 
S^^{q,uj) as given in eq. (j^) simplifies to 

Note that the 2;-component of the magnetization is conserved during the time evo- 
lution. Thus the neutron scattering function 5'|'^(q, cj) can be regrouped in terms of 
symmetry into the out-of-plane component S^^ and the in-plane component 

5f(q,a;) = 5f(q,a;), (11) 



with different physical behavior. (As reported in sections ^]3| and [4.4| , we find that for the 
two different components the exponent z is the same, but the scaling functions Q and / 
differ.) 



2.2 Dynamic Finite Size Scaling 

At the critical temperature Tkt and below, the XF-model is expected to be critical, 
with infinite correlation length C,- In this region the dynamic critical exponent z can be 
extracted by using the dynamic finite size scaling theory developed by Chen and Landau 

i- 

These authors also introduced a frequency resolution function in order to smoothen 
the effects of finite length of time integration, which was not necessary for most of the 
analyses in our study because of much longer integration times (see section ^). Their 
dynamic finite size scaling relations can then be simplified to 

^i^, = G'\u;L^'\gL) (12) 
L^5f (q) 

(replacing a factor u in front of S^{q,uj)), and 

c^^q, L) ^ u'J:{qL) = L-^"' n'\qL) , (13) 

analogous to eq. (|^) and eq. (pil|). 

We see that two different ways emerge to test dynamic scaling and to estimate the 
dynamic critical exponent z : Firstly, from eq. (0), z is given by the slope of a graph 
of logUm versus logL at fixed value of qL. Secondly, eq. (|1^) implies that for correctly 
chosen z and at fixed value of qL, graphs of (S'^'^(q, a;) / S'l''{q)) versus uL^ should all 
fall onto the same curve for different lattice sizes. Both procedures will only be valid for 
sufficiently large lattice size. 



3 



2.3 Analytical Results 



The dynamics of two-dimensional ferromagnets with easy-plane asymmetry, specifically 
the XY-mode\, were first analyzed by Villain [|^ and by Moussa and Villain [0. The 
in-plane scattering function was found to have a delta-function spin-wave peak at low 
temperature, and a spin-wave peak of the form 

^'■'■(g,^) ~ I (14) 



close to Tkt- Here rj{T) is the critical exponent describing the decay of the static spin-spin 
correlation function (and we now expect ?7 = | at Tkt)- 

Nelson and Fisher treated the classical XY model in a fixed length hydrodynamic 
description for T < Tkt, without vortex contributions. They obtained the transverse 
spin-spin correlation function, 

C-i^r,^-^,%{f), (15) 

for r, ct ^ 1 and ct ^ where c is the spin- wave velocity. The fourier transform of eq. 
(p!5|) has the form 

S^^iq^u) ~ ^vl/f^) , (16) 
q-^ ^ \cq J 

where the scaling function ^' behaves like 

1 1 ~ 1/ I 

around the spin-wave peak, and 

S^^{q,iu) r^Lu"-'^ (18) 

for large values of uo/cq. 

Nelson and Fisher also predicted that the dynamic critical exponent z, eq. (|^), which 
is expected to be 2; = d/2 for c? > 2, is 

z = I for < 2 . (19) 

Note that the value z = 1 and a linear dispersion relation are also implicit in the argument 
u/cq of the scaling function in eq. (0). 

Finally, both Villain and Nelson and Fisher predicted a very narrow (delta-function) 
spin- wave peak in the out-of-plane function S^^, at u = cq. 

More recently, Menezes, Pires, and Gouvea have performed a low temperature 
calculation which includes the contribution of out-of-plane fluctuations to the in-plane 
correlation functions. They worked in the harmonic spin-wave approximation which is 
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justifiable for large spin s and therefore also for our classical spins, and they used a 
projection operator technique. They found a spin-wave peak similar to that of Nelson 
and Fisher, 

nxx/ \ 2 1 T f e"^/** <^ k <^ TT and 

^ ~ " ?M|T^ • i e-V" « |1 - i^l « . . P") 

where u = uj/{cq). In addition to the spin- wave peak, they found a logarithmically 
diverging central peak, i.e. a signal at very small u, which diverges like 

S^^(q,uj) ~ — - + (less divergent terms) . (21) 

q log \uj\ 

Of course, a central peak at small temperature can also be caused by other mechanisms, 
for example vortex pairs diffusing like a dilute pair of solitons ||T0| . 

The dynamic behavior of the Xy-model is different above Tkt- For a phase transition 
of Kosterlitz-Thouless type, the spin stiffness should drop discontinuously to zero at Tkt, 
i.e. the spin- wave peak is predicted to disappear Above Tkt-, vortex-antivortex 

pairs unbind, and their diffusion leads to a strong central peak in S{q,uj). 

Mertens et al. [|l^, calculated S'(g, u) above Tkti assuming an ideal dilute gas of 
unbound vortices moving in the presence of renormalized spin-waves, and screened by the 
remaining vortex-antivortex pairs. They found a Lorentzian central peak for S^^, 

5""(g,cu) 2, (22) 

and a Gaussian central peak for S*^^, 

5-(g,^) ^ !^e-(^)' , (23) 

where 7 = ^y/nu/^,, it is the rms vortex velocity, ^ the correlation length, and n„ ~ (2^)~^ 
the free vortex density; and they compared their results to numerical simulations (see 
below). 



2.4 Previous numerical work 



Gerling and Landau [|ri|] carried out spin dynamics simulations on the XY model with 
L < 204 and found both spin-wave peaks and a central peak. The resolution was too 
limited, however, to allow quantitative comparison with theory or to extract an estimate 
for the dynamic exponent. 

Mertens et al. [0, [16[ performed spin dynamics simulations, with fixed system size 
100 and very low statistics (3 starting configurations) at T = 0.5 and T = 1.1. Below 



L 



Tc they observed only spin- wave peaks in both S^^ and S*^^; above Tc they saw a strong 
central peak and a weak spin- wave peak in S^^, and vice versa in S^^. The width and 



intensity of the central peaks were compatible with eq. (p2|) and eq 

Other earlier numerical work on dynamical behavior has been exploratory ||17| , |1 
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2.5 Experiments 



The closest physical realizations of the XY-model are materials with very large anisotropy, 
more generally described by strongly anisotropic Heisenberg models. Several experiments 
have studied the dynamics of such materials like Rb2CrCl4, K2CUF4, and CoCl2- 



In a recent study on stage-2 C0CI2 intercalated graphite, Wiesler et al. |0 found 
four temperature regimes with different behavior. There are indications of a Kosterlitz- 
Thouless transition at a temperature "T^j", though some properties disagree with KT 
predictions. Between temperatures "T/" and "T„", they observed spin- wave peaks. It 
is not clear whether ta central peak is present there. (In this temperature region the 
long range part of the scattering function shows true 2-dimensional character, whereas for 
T < Ti 3-dimensional correlations develop.) Above "T^", the in-plane scattering function 
showed the expected central peak, and the out-of-plane function exhibited damped spin- 
waves. 

In experimentally available materials both defects and the effects of residual three- 
dimensional couplings limit the effective size of the two-dimensional KT-like system to 
Lexp = 0(100) lattice spacings [0, |l9l. Remarkably, this size is similar to the lattice sizes 
of the present numerical study. For further discussion and an extensive listing of relevant 



literature, see the recent overview contained in 19 



3 Simulations 

We have studied the two-dimensional classical Xy-model with Hamiltonian given in eq. 
(0) on LxL lattices with periodic boundary conditions for 16 < L < 192, at temperatures 
T = 0.4, 0.6, 0.7, 0.725, and 0.8 in units of J/ks- Most of these temperatures are in the 
critical region T<Tkt = 0.700(5) (see appendix). 

Equilibrium configurations were created at each temperature using a Monte Carlo 
method which combined cluster updates of the x and y spin components (using the Wolff 
embedding method |^, ^) with vectorized Metropolis and overrelaxation spin re- 
orientations. After each single-cluster update, two Metropolis and eight overrelaxation 
sweeps were performed. Use of the cluster algorithm was important, since critical slowing 
down was severe for most of our simulations; the inclusion of cluster flipping reduced 
Monte Carlo autocorrelation times at L = 192 and T = 0.6 from more than 300 to about 
3 hybrid sweeps, while requiring only a factor of two more CPU time per sweep. We 
performed 200 hybrid sweeps between equilibrium configurations, and discarded the first 
5000 hybrid sweeps for equilibration. 

Between 500 and 1200 equilibrium configurations were generated for each lattice size 
and temperature. We found this many configurations to be necessary in order to suffi- 
ciently reduce statistical errors in the resulting neutron scattering function. The error 
bars in our figures represent the statistical errors for averages over the equilibrium con- 
figurations, drawn from the canonical ensemble. 
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Starting with each equihbrium configuration, the time dependence of the spins was 
determined from the coupled set of equations of motion, eq. (0), and was integrated 
numerically using a vectorized fourth order predictor-corrector method 0, with a time 
step size of 5t = O.OIJ^^. The maximum integration time was generally tmax = 400J~^; a 
few runs were also performed for lattice size 256 x 256 with tmax = 800 and produced 
the same physical results. 

The time-displaced, space-displaced spin-spin correlation functions C(r — r',t — t'), eq. 
(1), were measured for each time integration, with 

< t' < 0.1 tmax and < (t-t') < t,utoff = 0.9 tmax , (24) 

and were then averaged. By fourier transformation in space and in time, eq. (^, we 
obtained the neutron scattering function S{q,uj). The time integration in eq. (^) was 
performed using Simpson's rule, with a time step of O.IJ, which has been shown to be 
sufficiently small. 

To reduce memory and computer time requirements, we restricted ourselves to mo- 
menta q = (g, 0) and (0, q), with q determined by the periodic boundary conditions, 

27r , , 

q = ng—, Uq = 1,2,..., L, (25) 

and data from these two spatially equivalent directions were averaged together to further 
enhance the statistical accuracy. We used fast fourier transforms to increase the efficiency 
of the program in calculating correlation functions. 

The frequency resolution Auj of our results is determined by the time integration 
cutoff tcutoff = 0.9 tmax, sce eq. (0), which will introduce oscillations of period 2n/tcutoff 
into S{q,uj). Since we observed very sharp spin-wave peaks (see section [4.1|) , we chose 
to integrate the equations of motion to very large times. We used tmax = 200J^^ for 
L < 96, and tmax = 400J^^ for L > 128. A theoretical delta-function in frequency will 
then become a widened peak with a width at half maximum of 

in the simulation data. To smoothen the oscillations, previous spin dynamics studies 
ni , H have employed a frequency resolution function, replacing 



C^^(r,t) by C^^(r,t) exp{-^{t Scof) (27) 

to compute S{q,uj). Because of the large values of tmax in our study, we achieved a 
very small frequency resolution Aa;, and the oscillations were not noticeable for most of 
our data. We therefore did not generally use a frequency resolution function and could 
significantly simplify the analysis of our data. 
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The numerical integration of the equations of motion can potentially become unsta- 
ble at very large integration times. We checked that for our calculations, in which we 
integrate to much larger times than previous studies, we do not encounter this problem. 
We verified that the constants of motion (energy and magnetization in z-direction) do 
remain constant, with a relative variation of less than 3 x 10~^. We also verified that the 
neutron scattering function remains virtually unchanged when an additional integration 
of length t = 200J~^ is performed from each equilibrium configuration before starting 
to calculate time-displaced spin-spin correlation functions. All simulations were carried 
out using highly vectorized programs on the Cray C90 at the Pittsburgh Supercomputing 
Center. 

4 Results 

We now present our results for the dynamic structure factor S{q,u), its dependence on 
temperature, frequency, momentum, and lattice size, and we analyze its dynamic scaling 
behavior. With few exceptions we have analyzed the data without the use of a frequency 
resolution function eq. ([27|). (The effects of such a function, and of integrating to shorter 
maximum times are described together with initial results in reference ||25|| .) 

In order to investigate the XF-model in the critical phase, we chose several tem- 
peratures T < Tkt, including the best previous estimate of T^t ~ 0.725J/kB 0, and 
one temperature well above the transition, T = 0.8 J/ ks- The results of our analysis 
prompted us to perform additional static Monte Carlo studies, which are described in the 
appendix. They provided an improved estimate of T^t = 0.700(5) J/fc^. In order to elu- 
cidate the situation closer to the transition, we performed additional (but less extensive) 
spin dynamics simulations at T = 0.700 J/ ks- 

4.1 Spin- Wave Peak 

Figure 1 shows the temperature dependence of S{q,uj) as a function of u, for lattice 
size L = 128 and fixed small momentum q = tt/A8 (i.e. ra^ = 2 in eq. (p5|)). Here, as 
in other results which we shall show, the error bars are determined from the statistical 
variation of results obtained from different initial spin states. The in-plane component 
gxx^ fig. 1(a), exhibits a very strong and moderately sharp spin- wave peak at temperatures 
T^Tkt- Even at the lowest temperature, however, the width of the peak is larger than the 
minimum value eq. (|26|) due to finite cutoff time. The position of the peak moves towards 
lower uj as the temperature increases, and the peak broadens slightly. Just above the 
transition, at T = 0.725, there is still both a strong spin- wave peak and a sizeable central 
peak in S^^. At higher temperature, the spin- wave peak disappears completely (for this 
low momentum) and only a large central peak remains. Note that from KT-theory |^ 
one would expect complete disappearance of a spin- wave peak at all T>Tkt- 

There is sizeable additional structure in 5*^"^ away from the spin-wave peak at temper- 
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Fig. 1. Temperature dependence of the neutron scattering function 
S{q,uj) as a function of frequency uj. The transition temperature is 
Tkt ~ 0.700(5). Lattice size L = 128 and momentum q = 2 x ^ = 
in all cases, (a) xx- component, (b) 2;z-component; the inset shows the 
data at T = 0.6, smoothened with a resolution function, eq. (pTf), with 
Slu = 0.015. 



9 



atures up to Tkt- We will discuss this structure in the following sub-section. 

The out-of-plane component S*^^, shown in fig. 1(b), is two orders of magnitude weaker 
than the in-plane component. It exhibits a very sharp spin- wave peak for T < Tkt-, 
whose width is limited by our tu-resolution. The finite time cutoff eq. (^) produces very 
noticeable oscillations in S'(g, uo). (The magnitude of these oscillations is minute compared 
to the intensity of the spin- wave peaks in S^^ .) The oscillations can be smoothened by 
convoluting S{q,uj) with a gaussian resolution function in frequency, as is shown in the 
inset. No central peaks are visible in S^^ at T < Tkt- At T = 0.725, the peak in S^^ 
is still present, with a larger width similar to that in S^^. In contrast to 5*^"^, there is a 
clear, but weak, spin-wave peak in S^^ even at T = 0.8 > Tkt and small momentum. It 
is of similar intensity as the peak at lower temperatures. (See also section 2.5). 

Figure 2 shows the lattice size dependence of S{q,uj), at fixed momentum q = n/lG. 
Below the transition. Fig. 2(a), the intensity of the spin- wave peak depends strongly 
on lattice size, whereas its position is constant. (The out-of-plane component S^^ is 
dominated by finite time cutoff effects for T < Tkt, and we do not show it here). Just 
above the transition, at T = 0.725, the spin-wave peak in 5*^^ appears to gain intensity 
slightly as L increases, whereas neither its central peak nor the spin-wave peak in S^^ 
show any finite size effects. At higher temperature, fig. 2(d) and 2(e), there is no visible 
lattice size dependence in either S^^ or S^^. Notice the two different vertical scales for 
S^^ and S^^. Data taken for L = 16 and L = 32 exhibit such strong finite size rounding 
that we have chosen not to show the data here. 

In figure 3 we show the momentum dependence of S{q,uj). Fig. 3(a) and (b) display 
the behavior at T = 0.725, which is qualitatively similar to that at lower temperatures. 
The position of the spin-wave peak is the same for 5*^^ and S^^ and is proportional to 
momentum for small q. As q increases, the peak broadens, and becomes less intense, yet it 
remains well defined. The additional structure in S^^ is strongly momentum-dependent, 
as will be discussed below. 

For the 2;2;-component, both the total intensity and the relative loss of intensity with 
increasing momentum are much smaller. Our a;-resolution dominates the width of the 
spin-wave peak in S^^ only at the smallest q (which also appears in figure 1); it is not 
dominant at higher momenta or for S*^^. We conclude that S^^ has the expected delta- 
function form only for very small momentum. (Higher order perturbation theory also 
predicts a finite width |2^.) As shown in fig. 3(c), the intensity of the spin- wave peak 
decreases much more rapidly for S^^ than for S^^ with increasing g; the intensities cross 
each other well before the zone edge is reached. This behavior is similar at other tem- 
peratures. We also note that at all temperatures the total intensity S^^{q) (not shown) is 
constant with q, whereas S^^{q) decreases (as g^"^) and crosses S^^{q) at a slightly larger 
momentum q ^ 1.5 — 2. 

Well above the transition, at T = 0.8, fig. 3(d), S^^ has no noticeable spin-wave 
peak at small momentum. The strong central peak rapidly loses intensity with increasing 
momentum. In marked contrast, the behavior of S^^ (fig. 3(e)) is very similar to that 
at lower temperatures, with clear but broadened spin-wave peaks. Notice that there is 
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(S3/J 

Fig. 2. Lattice size dependence of S{q,u!), at fixed momentum q — 7r/16. 
(a)-(d): xx-component; (e): 2;2;-component 
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Fig. 3. Momentum dependence of S{q,uj) at fixed lattice size L — 192. 
(a) T=0.725, xa^-component; (b) T=0.725, zz-component; (c) T=0.7, 
height of spin- wave peak as function of g (at L = 128); (d) T=0.8, xx- 
component; (e) T=0.8, 2;2;-component; (f) T=0.8, xx-component at large 
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now non-zero intensity at small a; in S^^. Remarkably, at very large momenta spin- waves 
appear in S^^ even for T = 0.8 (fig. 3(f), so that both a central peak and a spin-wave 
peak are present. Note that the vertical scale in fig. 3(f) is 100 times smaller than in fig. 
3(d). There is no noticeable lattice size dependence here. 




Fig. 4. Dispersion relation: spin-wave frequency against momentum, 
at L = 192 and four different temperatures. Note that for T = 0.8 only 
S^^ has a noticeable spin- wave peak, S^^ does not. The inset shows the 
spin-wave velocity |p as a function of temperature. 



Figure 4 shows the position ujp of the spin-wave peak as a function of momentum. 
The expected linear portion of the dispersion curve extends to rather large momenta. 
With increasing temperature, the spin- wave velocity dujp/dq, which is proportional to the 
spin- wave stiffness, decreases slowly and approximately linearly, as shown in the inset. 



and theoretically expected for small T [p!4| , At T < Tkt, i^p is the same for S^^ and 
for S^^, as expected by theory 0, |[. At T = 0.8 on the other hand, we can only plot 
the position of the residual peak in 5*^^, because 5*^^ has dropped sharply to zero here, as 
expected for a KT transition. 
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4.2 Additional Structure in S^^{q,u) 



If we expand the vertical scale in plots of S^^{q,uj), we find that the in-plane component 
S^^{q, u) shows rich structure in addition to the spin- wave peak. Note that the intensity of 
this structure is typically 10~^ of the maximum. It is visible at all temperatures T^Tkt- 
At the lowest temperature, the absolute intensity of this structure is low, but the relative 
intensity is quite high (see also fig. 11(a)). At higher temperature, the structure becomes 
rather smeared. No such structure can be found in S^^{q,uj). The locations of 

additional maxima in S^^ are essentially unchanged with L when = gL/(27r) is held 
fixed. Figures 5(b) and (c) show, on a logarithmic scale, that for odd values of there 
are strikingly regular pronounced peaks at = ojp/uq, 2,(jjp/nq, and for even such 
peaks appear at = 0, loOpluq, Aup/nq, .... At large ra^, fig. 5(d), individual peaks cannot 
be distinguished; instead S^^ is nearly constant below the spin-wave peak there. (In figs. 
5(b) and (c), the data have been smoothened slightly with a resolution function eq. (|27|) 
with 5uj = 0.01, in order to reduce the wiggles and allow general features to be identified.) 
In addition to the regularly spaced pronounced peaks, there is further "fine structure" in 
S^^, clearly visible in fig. 5(a). In the course of our study, the additional structure became 
clearer as the statistical quality of the data improved; it is apparent that the structure is 
not statistical noise. Very close to the spin-wave peak, part of the additional structure 
may be due to the finite-time cutoff in our time-integration; but most of the observed 
structure must be due to different reasons. 

One simple explanation for the observed rich structure, which is consistent with the 
data but for which we have no rigorous theory, is that of multi-spin-wave effects. Of 
these, two-spin-wave processes are likely to be the most important. Thus, at a given total 
momentum q we can have either a single spin-wave excitation of momentum q, or two 
spin-waves for which the sum or difference of momenta equals q. The result will then be 
both a single-spin- wave peak at a characteristic frequency ujp{q) as well as additional sum 
and difference peaks due to the two-spin-wave processes. 

Of particular interest is the case when the two spin-waves have momenta qi and q2 
that are coUinear, so that g = gi -|- g2 is a scalar equation. Since the momenta are discrete 
on a finite lattice, qi = Ug.^, this implies Uq = Uq^ +nq^. With a linear dispersion relation 
LJ = cq, the difference of the two spin-wave frequencies is then co = {2nq^ — Uq) c^, i.e. 
just the series of additional peaks that are visible in figures 5(b) and (c). 

Using the measured dispersion relation (fig. 4), we have calculated the frequencies of 
two-spin-wave excitations consisting of the most likely individual spin-waves, i.e. those 
with smallest individual momenta. For the case of = 3 and T = 0.6, these locations 
are marked in fig. 5 (e) . They are identified by the coordinates Uq-^ of one of the spin- waves 
in reciprocal space; the sum of the two spin-wave momenta must equal q= (3 x ^,0). 
The locations of the resultant excitations agree extremely well with the positions of the 
small peaks in S{q,uj), but we have no way of comparing intensities. 

The presence of distinct small peaks at = 0, 2ujp/nq, iup/uq, .... for even Uq at 
T < Tkt complicates the identification of a possible central peak. Interpolating the 
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Fig. 5. High resolution study of the "fine structure" in S^^{q,uj) for 
T = 0.6, L = 192. (a) Ug = 4; note that the maximum value of the spin- 
wave peak is ~ 160. (b),(c) log plot of S^^{q,u) at small values of n^; 
the data are smoothened, with 6u; — 0.01; (d) S^^{q,uj) at q — 7r/2; (e) 
Ug = 3: vertical arrows show the location of two-spin-wave peaks formed 
by spin- waves of small momentum |^ < 4 
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intensities for odd values of Uq (which do not show peaks at = 0) to obtain estimates 
for even Uq, we conclude that there is indeed extra intensity at a; = which is not 
attributable to two-spin-wave processes. 

4.3 Finite size scaling of characteristic frequency 

Equation (||) defines the characteristic frequency Um of the whole spectrum of S{q,u!). 
When there is only a single spin-wave peak, then Um coincides with the spin-wave fre- 
quency Up. This is the case at T = 0.4, where all frequencies coincide (within error bars), 
= = = uj^. Closer to the transition, intensity between = and the 
spin-wave peak grows; therefore the characteristic frequency uj^ becomes smaller than 
the spin-wave frequency Up. Their difference is thus a measure of the relative weight of 
non-single-spin-wave contributions to S{q,u!). Figure 6(a) shows the situation at Tkt, 
where 5*^^ exhibits large non-single- spin- wave contributions and oo^ < ujp^. 




Fig. 6. Characteristic frequency ujm and spin- wave frequency ujp for S^^ 
and S"^^, at L = 128. (At T = 0.8 there is no xx-spin-wave frequency 



Above the transition, a;^^ is no longer linear in momentum for small q, as shown in fig. 
6(b), and differs strongly from the 2;2;-component. The latter still has both a spin-wave 
peak that is linear in momentum, as well as intensity at small lj, so that is smaller 
than Up^. 

The central question of critical dynamics is that of scaling, i.e. whether data from 
lattices of different size match when properly scaled. As mentioned in section ( p^.2| ), we 
can test scaling and extract the dynamic critical exponent z in two ways, by analyzing 
the characteristic frequency Um, or by looking at S{q,uj) itself. 

We concentrate on the dynamic critical behavior of our model at T < Tkt, the critical 
region, in which the correlation length in an infinite system is divergent. The relevant 
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length scale on a finite lattice is therefore the lattice size L, and we expect scaling for 
suitable functions of qL, as described in section ( p.2|) . From the analytical results we 
expect the dynamic critical exponent to be 2; = 1 (eq. ([T9|) ). 

In figure 7, (a),(b) we show uj^ U as a function of qL, for z = 1.00 and at temperatures 
T < Tkt- From eq. (|13D we expect the data to fall on a single curve if scaling holds. This 
is indeed exactly the observed behavior at all temperatures T < T^t- The asymptotic 
behavior for large L is strictly linear, UmL^ ~ qL; i.e. for z = 1, Um ~ q- For 




Fig. 7. Finite size scaling of the characteristic frequency; is plotted 

against gL, for z = 1.00. 

each finite lattice size the dispersion curve flattens when q becomes large. Therefore as 
L increases, the data in figure 7 start to move away from the asymptotic behavior at 
progressively larger values of qL. We have also analyzed the scaling behavior for different 
values of z, e.g. z = 1.10 (not shown in the figures). In that case the data for different 
L diverge from each other immediately; they do not fall onto a common line even at the 
smallest momenta. 

The scaling curves for at all three temperatures T < T^t are very similar, with 
variation only in their slope. In contrast, we do not observe similar scaling behavior in 
at T = 0.8 above the transition (not shown in the figures). 

Analyzing the out-of-plane characteristic frequency u^, we found that (for g 7^ 0) it 
has the same scaling behavior as the in-plane component. At T = 0.4 the data for uj^ 
and are indistinguishable. When intensity below the spin-wave peak grows in S^^ at 
larger T, the scaling curve for uj^ has a smaller slope than u^. Interestingly, at T = 0.8, 
above the transition, not only are there spin- wave peaks present in S^^, but uj^ also shows 
the same scaling behavior as below the transition, with z = 1.0. 
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4.4 Finite size scaling of S{q,u) 



If dynamic finite size scaling holds, then the scaled neutron scattering function itself 



should fall onto a single curve for sufficiently large lattices. Corresponding to eq. (p!2D, 
figure 8 shows S^^{q,uj)/ {L^S{q)) versus uL^. 

We see that at all temperatures T^Tkt the data do indeed fall onto a single line 
within error bars, when we choose z = 1.00. This is not only true for the spin-wave 
peak itself, but for the whole range of u L^. Only the data from very small lattices (not 
shown here) deviate systematically. Even at T = 0.725 the data scale quite well for the 
values of L for which we have data. The correlation length at T = 0.725 is still very large 
(see appendix); deviations from scaling could presumably be seen if data for much larger 
lattices were available. 

Note that scaling with u implies that at fixed qL and for large lattices the spin- 
wave peak is very narrow in units of u. Its width is therefore very sensitive to the time 
cutoff in the spin dynamics integration, and we had to use the very long time integrations 
described in section in order to obtain asymptotic results. 

The finite size scaling behavior is very sensitive to variations in z. As an example, 
figure 8(e) shows that at T = 0.4 the data do not scale when choosing z = 1.05, even 
upon excluding all lattice sizes L < 96. Using similar plots, we obtain 



(2^ 



(The relatively large error for T = 0.7 is a consequence of the limited amount of data 
available at this temperature.) It is remarkable that the dynamic critical exponent is the 
same accross this range of temperatures, whereas the static exponent rj varies strongly, 
from T] = 0.082(2) at T = 0.4 to r] = 0.247(6) at T = 0.7 (see appendix). 

The z;z-component of S{q, u) is extremely narrow at T = 0.4 and T = 0.6, and cannot 
show scaling given our maximum integration times. At T = 0.7, the spin- wave peak in 
S^^ has become wider, and we do observe scaling, as shown in figure 8(f). 

Above the phase transition, the relevant length scale is the correlation length ^, not 
the lattice size L. We saw earlier (fig. 2(d),(e)) that finite size effects are already small for 
our lattice sizes. Yet, surprisingly, there is scaling-like behavior for small momenta in S^^ 
even above the transition, as show in fig. 9(a). Note that, at constant n^, the horizontal 
scale uL is proportional to u/q. The data do not scale when different Ug are compared. 
At large momenta, a spin- wave peak is visible (see also fig. 3(f)). For the out-of-plane 
component 5"^^ (fig. 9(b)), the spin- wave heights do not obey the scaling equation, eq. 
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Fig. 8. Finite size scaling of the neutron scattering function; 

S{q,Lj)/{L^S{q)) is shown versus ojL^, with constant nq = ^ = 2. The 
critical exponent is 2; = 1.00 in (a)-(d) and (f). (a)-(d): 5"^^, z — 1.00; 
(e) 2; = 1.05 for comparison; (f) z^^-component. 
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Fig. 9. Finite size scaling plot of S{q,uj) for T = 0.8, above T^t, with 
constant rig = ^ =2. (a) 5^^, (b) 



4.5 Tests of Nelson-Fisher scaling form; large cj-behavior 

Nelson and Fisher [|| predicted the scaling form eq. ( |T^ ) for S^^{q,uj) : 




cq 



This provides an explicit opportunity to compare data at different temperatures and at 
different values of Ug. Equation ([T^) implies the finite size scaling equation ([T^) used in 
the previous section, with z = 1: 



S''%q,u) 



\cq 



LS'^'^iq) cqL f\iJ (^] 

■I \cqj cq 



f{qL, uL) . (29) 



It also implies 




for large enough lattice sizes L, for which eq. (|TB|) is valid. (Note that when rig = |^ is 
constant, the arguments uL in eq. ( p9D and ^ in eq. (^) are equivalent.) 

In fig. 10(a) we use eq. (|30|) to compare the in-plane scattering function S'^^ at different 
temperatures, for constant Ug = 2. The same data appear unsealed in fig. 1(a), and with 
dynamic finite size scaling in fig. 8. Obviously, eq. (^0|) is not satisfied: the scattering 



function at different temperatures within the KT-phase does not scale to the same shape 
This is also the case at other values of n, 



-cq' 



q- 



Different values of Ug are compared in fig. 10(b), at T = 0.7. Again, the data do not 
scale. Moreover, the dependence of the spin-wave peak on Ug is not monotonous: for 
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Fig. 10. Tests of the Nelson-Fisher scahng form, eq. ([Tq) . (a) Different 
temperatures; (b) Different values of n^. 



increasing the peak height first grows, is approximately constant for = 3... 8, and 
then shrinks. As becomes large, there is growing intensity below the spin-wave peak 
(see also fig. 5(d)). Note that equations (p!6| ) and (|30| ) are at odds with the fact that 
the additional peaks we observed in 5*^^ (section |4.2|) have positions which do depend on 
nq. The data in fig. 10(b) have been obtained with constant time-cutoff tcutojf = 360. 
The picture is virtually unchanged when data with t cutoff ~ ^ (and t cutoff sufficiently 
large) are used. A similar comparison at other temperatures T < T^t shows still stronger 
deviations from scaling. 
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Nelson-Fisher scaling form, eq. (|16]). The data are smoothened, with 
6u = 0.015. (a) Different temperatures, (b) Different values of Ug. 
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For large frequencies, S{q, u) appears to be independent of lattice size. We show S'^^on 
a log-log scale in fig. 11, with data scaled similarly as in fig. 10. Nelson and Fisher |^ 
predict that S^^ is also independent of momentum and follows a power law uj~^, with 
p = 3 — f] (eq. (0)). We see somewhat different behavior, which is determined mainly 
by Ug. For rig = 2, fig. 11(a), the data can be fitted with p = 3.0(1) at all temperatures. 
Note the sizeable structure at low T. For larger Ug, fig. 11(b), p increases. There is 
also noticeable curvature in u, with larger p at higher uj; p also increases slightly with 
temperature. The exponents in fig. 11(b) range from 3.7(1) at Ug = A to 5.4(2) at 
Ug = 24. The out-of-plane correlations S^^ (not shown) also decay with a power law with 
momentum-dependent exponents. 

Equation (|1^) also implies 

q'-^S'%q,u) = ^(^) • (31) 

In fig. 12 we use eq. (|3TD to compare data for different momenta q, at the KT phase 
transition temperature, using r] = 0.25 and constant rig. Here the data do scale. This 
scaling is also implied by dynamic finite size scaling, fig. 8(c), together with a functional 
dependence S^^{q) ~ g''"^. Note that for constant L (instead of constant rig), we obtain 
a non-scaling picture indistinguishable from fig. 10(b). 
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Fig. 12. Test of the Nelson-Fisher scaling form, for different lattice sizes 
L. 



4.6 Lineshapes 

In figure 13 we compare our results with theoretical predictions for the shape of S^^{q, u). 
We show data at T = 0.7, for L = 128 and g = 7r/32 , normalized according to eq. (PP]), 
and we compare with predictions using 77 = 1/4, similarly normalized. 
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Fig. 13. Comparison of the lineshape of S^^{q,uj) with theoretical pre- 
dictions. Data are at T = Tkti L = 128 and q = 7r/32 (thick hne), 
and are normahzed according to eq. (|30D. The two thin hnes represent 
the predictions by Nelson and Fisher eq. (|T6|) (continuous line) and by 
Villain eq. (p^ ) (dashed line), both with rj = 0.25. The thick dashed line 
is obtained by a fit to the data with an ad-hoc function related to the 
Nelson-Fisher form (see text). The inset shows the data and predictions 
on a log-log plot that includes large values of u. 



The predictions by Nelson and Fisher and by Villain both have a pole at the spin- 
wave peak (eqs. and |1^, as shown in the figure. In order to compare better with our 
data which have been obtained with a time integration of finite length, we also tried to 
convolute the predictions with the fourier-transform of a finite-time cutoff at t^utoff = 360. 
The resulting functions (not shown) exhibit very strong oscillations (size= 0(0.5)), and a 
spin- wave peak that is much higher (about 3.7) and more narrow than that of the data. 
We conclude that S^^{q,t) decays much faster in time than predicted. 

As mentioned before, the rich structure in S^^{q,uj) below and above the spin-wave 
peak that was described in section had not been predicted, except for a (small) central 
peak, eq. (|2T|) , foreseen by Menezes et al. [^. 

Since the Nelson-Fisher prediction does not agree well with the data at the spin-wave 
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peak, we tried to find a functional form which does fit these data reasonably well and 
should thus be an approximation to the actual form, so that our data can more easily 
be compared to the results of future theoretical calculations. We found that a modified 
form of eq. ( ]30| ) works well, namely a widening of eq. (|30| ) with a Gaussian resolution 
function eq. (|27|) , with 5uj a free parameter. Around the spin- wave peak we obtained 
fairly good agreement with our (unconvoluted) data, as shown by the thick dashed line 
in Fig. 13, which uses 5uj = 0.01. However, different values of 6uj in the modified function 
are necessary to describe the data at different n^. The modified function can of course 
not describe the additional structure in S^^{q,uj), including the central peak. 

The large frequency behavior of the data is shown in the inset in figure 13. The 
prediction by Nelson and Fisher agrees with the data at large frequencies qualitatively. 
However, as described in section ^]5|, a fit to S^^{q, u) for large u results in a different power 
law exponent than predicted. Not surprisingly, at large frequencies the prediction by 
Villain, intended for the spin-wave peak divergence, does not describe the data correctly. 

We conclude that below and above the transition, the actual lineshape is quite different 
from the predicted forms, with a much wider spin-wave peak, and a lot of additional 
structure. 

Above the transition, the theoretical predictions eq. (^) and eq. ( p3D do not describe 
the data well either. For the in-plane component S^^ we see two different regimes in a;. At 
small uj, it is compatible with a Lorentzian-like peak ~ {u^ + c)"^, but with an exponent 
b that is momentum-dependent (e.g. b ^ 1-1(1) at g = 7r/48, b ~ 0.43(2) at g = 7r/16). 
At large u, S^^ decays with a power law ~ cj"^, with c = 3.2(2) for small momenta (see 
fig. 11(a)). The out-of-plane component S*^^ does not show the predicted central peak at 
all; instead it exhibits a spin-wave peak. 

5 Conclusions 

We have performed the first high precision study of the dynamic critical behavior of the 
XF-model, at five different temperatures below, at, and above Tkt^ on square lattices of 
size up to 192 X 192. We have determined the critical temperature to be T = 0.700(5) J/ks- 
Starting from about 1000 equilibrium configurations generated by an efficient Monte Carlo 
procedure at each temperature and lattice size, we have integrated the equations of motion 
of the spins to very large times, tmax = 400 and measured space-displaced, time- 
displaced correlation functions to compute the neutron scattering function S{q,uj). 

At temperatures up to T^t, S{q, u) exhibits very strong and sharp spin-wave peaks 
in the in-plane-component S^^. As T increases, they widen slightly and move to lower 
UJ, but remain pronounced even just above T^t- For increasing momentum they broaden 
and rapidly lose intensity. Well above Tkt, the spin-wave peak disappears in S^^, as 
expected, and we observe a large central peak instead. 

In addition to the spin-wave peak, the in-plane component 5*^^ exhibits a rich structure 
of small intensity, which we interpret to come from two-spin-wave processes. Furthermore, 
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gxx sJiows a clear central peak, even below Tkt-, which becomes very pronounced towards 
the critical temperature. Neither this strong central peak nor the additional structure are 
predicted by existing analytical calculations. 

The out-of-plane component S^^ is much weaker than S^^, except for large q. It 
displays a sharp spin-wave peak at all temperatures, even above Tkt- The peak widens 
with increasing temperature, and only at low T is it consistent with a delta-function 
shape. 

We measure the dispersion relation, i.e. the position of the spin-wave peak as a function 
of momentum, to be linear to high accuracy. Its slope, the spin-wave velocity, decreases 
with increasing temperature approximately linearly, as expected from approximate ana- 
lytical results at small T. 

Examining dynamic finite size scaling, we show that both the characteristic frequency 
ujm and the neutron scattering function S'(g, uo) itself scale very well for all T < Txt, with 
a dynamic critical exponent oi z = 1.00(4) that does not depend on temperature, whereas 
the static exponent 77 varies strongly. 

The shape of the scaling function is not well described around the spin-wave peak by 
the available theoretical predictions, nor is the shape of the scattering function above the 
transition, and the additional structure had not been predicted at all. The data which 
we have presented here are of sufficiently high quality that meaningful comparison with 
theory and experiment is possible. We hope that this spin dynamics study will thus serve 
to stimulate further effort in this area. 
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Appendix: Static critical behavior 



The determination of the transition temperature in XF-hke systems has been notoriously 
difficult 0, 1^. The best previous estimate pf for Tkt for the model considered here 
was Tkt = 0.725 ± 0.010, estimated from the onset of vortex-pair unbinding, which is 
a procedure that is quite difficult to apply with high precision. The results of our high 
resolution spin dynamics study at T = 0.725 prompted us to perform a new, more accurate 
determination of Tkt, using the powerful hybrid Monte Carlo algorithm described in 
section ^. 

We carried out a set of static Monte Carlo simulations, with lattice sizes L = 64, 128, 
and 256, and 40000 hybrid Monte Carlo sweeps in each run. We then analyzed the static 
correlations C(r) = ^{S^{0)S^{r) + 5'^(0)5'^(r)) in three different ways: (i) using finite 
size scahng, (ii) with a ffi to a power law decay, and (iii) with a fit to the free lattice 
propagator. The results of all three methods are in excellent agreement. 

The finite size scaling ansatz for the correlation function is 



C{r,L) = r-^ f{iL'/\^), (32) 

where i is the reduced temperature (t = |1— T/Td), and u the correlation length exponent. 
Since our model is critical throughout the KT-phase, we have t = for all T < Tkt- With 
the correct value of rj, the data for different lattice sizes should therefore coalesce on a 
plot of C(r) L'' versus ^. Figure 14 shows such plots at T = 0.7, for r] = 0.24, 0.25, and 
0.26. Choosing rj = 0.25 results in very good scaling over the whole range of distances r. 
At lower temperatures, scaling of similar quality is achieved with smaller values of rj. At 
a slightly higher temperature of T = 0.71, however, we observed only mediocre scaling, 
with an effective r] of 0.29(1), but with small systematic deviations from scaling already 
visible. At T = 0.725, the deviations from scaling are stronger still. 

Independent estimates of t] were obtained from simple power law fits C{r,L) ~ r~^, 
for lattice size L = 256, and Vmax < 20. The results for rj depend only very little on the 
fit-range, and agree within error bars with those from finite size scaling. 

A simple power law ansatz ignores the fact that C (r, L) is actually periodic in r with 
period L. The full functional form can be deduced from the fact that within the KT- 
phase, the model is thought |[l|, |2^ to behave like a free field theory, for which the exact 
finite lattice propagator (correlation function) at an effective "temperature" l/rj is [p9| 



C(r,L) = exp(-r/r(r,L)), with (33) 
27r ^ l-cos(rg^^) 
^ - ^,.ir^o4-2cos(g.f)-2cos(g,f)- ^^^^ 

We also used this functional form to fit C (r, L) . The results for different lattice sizes and 
different fit ranges (up to r = L/2 and excluding r < 5) agreed with each other and with 
the results from finite size scaling. The quality of ffis was very good for all T < 0.70, 
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Fig. 14. Finite size scaling plot of the static correlation function C^^ir) 
at T = 0.7 ~ T/^T; using three different estimates for r/. 



whereas for larger temperatures (T = 0.71 and above) it deteriorated strongly, and the 
results became lattice size dependent. 

We have found no evidence for logarithmic corrections. It is possible that small cor- 
rections are present and introduce a subtle bias. In this case our results for 77 would need 
modification. Assuming that logarithmic corrections are indeed negligible, we obtain as 
our combined results from all three methods 

r/ = 0.082(2) at T = 0.4 

r/ = 0.153(5) at T = 0.6 (35) 
77 = 0.247(6) at T = 0.7 . 

We note that the expected linear dependence of on T below Tkt does not seem to be 
satisfied at these temperatures in our model. 

Both from assuming ri = 1/4 at |jl|, and from the different qualitative behavior of 
C(r, V) for T > 0.710, we conclude that the Kosterlitz Thouless transition temperature is 

Tk^ = 0.700(5) J/Zcb. (36) 
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This estimate is slightly below the value from |^ which we had used at the start of the 
spin dynamics study, and clearly below the estimate of ref. . 

Slightly above the KT-transition, at T = 0.725, the dynamic behavior of S{q,uj) in 
our study resembles that at T^t (whereas at T = 0.8 it is very different). This may be 
explained by looking at the correlation lengths: The correlation function at T = 0.725 
can be fitted by an Ornstein-Zernike form 

r(r) ~ e-'^/« (37) 

with a value of ^ = 0(400). Since this correlation length is larger than the lattice sizes 
we have used in our study, a behavior resembling the KT-phase (where ^ = oo) is not 
surprising. A similar (approximate) fit at T = 0.8, on the other hand, gives ^ = 0(10), 
much smaller than our lattice sizes. After our work was completed, we received a paper 



by Cuccoli et al. with a Monte Carlo study of the statics of our model. Their results 



are in excellent agreement with ours. 
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